arXiv:1501.04836vl [cs.SC] 20 Jan 2015 


Subtropical Real Root Finding 


Thomas Sturm 

Max Planck Institute for Informatics 
Saarbriicken, Germany 
sturm@mpi-inf.mpg.de 

January 20, 2015 


Abstract 

We describe a new incomplete but terminating method for real root 
finding for large multivariate polynomials. We take an abstract view of 
the polynomial as the set of exponent vectors associated with sign infor¬ 
mation on the coefficients. Then we employ linear programming to heuris- 
tically find roots. There is a specialized variant for roots with exclusively 
positive coordinates, which is of considerable interest for applications in 
chemistry and systems biology. An implementation of our method com¬ 
bining the computer algebra system Reduce with the linear programming 
solver Gurobi has been successfully applied to input data originating from 
established mathematical models used in these areas. We have solved sev¬ 
eral hundred problems with up to more than 800000 monomials in up to 
10 variables with degrees up to 12. Our method has failed due to its 
incompleteness in less than 8 percent of the cases. 


1 Introduction 

Our work discussed here is motivated by our studies of Hopf bifurcations usi 
[U for reaction systems in chemistry and gene regulatory networks in systems 
biology, which are originally given by systems of ordinary differential equations. 
Hopf bifurcations can be described algebraically [3 EH HU [TO], resulting in 
one very large multivariate polynomial equation / = 0 subject to few much 
simpler polynomial side conditions gi >0, ..., >0. For such systems one is 

interested in feasibility over the reals and, in the positive case, in at least one 
feasible point. It turns out that, generally, scientifically meaningful information 
can be obtained already by checking only the feasibility of / = 0, which is the 
focus of this article. For further details on the scientific background, we refer 
the reader to our publications EHIIMIIMIEIIH]. 

With one of our models, viz. Mitogen-activated protein kinase (MAPK), we 
obtain and solve polynomials of considerable size. Our currently largest instance 
mapkeSeS contains 863438 monomials in 10 variables. One of the variables 
occurs with degree 12, all other variables occur with degree 5. Such problem 
sizes are clearly beyond the scope of classical methods in symbolic computation. 
To give an impression, the size of an input file with mapkeSeG in infix notation is 
30 MB large. HT[5]X-formatted printing of mapke5e6 would fill more than 3000 
pages in this document. The MAPK model actually yields even larger instances, 


1 



which we, unfortunately, cannot generate at present, because in our toolchain 
Maple cannot produce polynomials larger than 32 MB. 

This article introduces an incomplete but terminating algorithm for finding 
real roots of large multivariate polynomials. The principle idea is to take an 
abstract view of the polynomial as the set of its exponent vectors supplemented 
with sign information on the corresponding coefficients. To that extent, out 
approach is quite similar to tropical algebraic geometry |30j . However, after our 
abstraction we do not consider tropical varieties but employ linear programming 
to determine certain suitable points in the Newton polytope, which somewhat 
resembles successful approaches to sum-of-square decompositions [2^. 

We have implemented our algorithm in Reduce m using direct function calls 
to the dynamic library of the LP solver Gurobi m- In practical computations 
on several hundred examples, our method has failed do to its incompleteness in 
less than 8 percent of the cases. The longest computation time observed was 
around 16 s. As mentioned above, the limiting factor at present is the technical 
generation of even larger input. 

In Section|^we introduce a specialization of our method that only finds roots 
with all positive coordinates. This is highly relevant in our context of reaction 
networks, where typically all variables are known to be positive. We also discuss 
an illustrating example in detail. Section generalizes our method to arbitrary 
roots. In Section]^ we discuss issues and share experiences related to a practical 
implementation of our method. In Sectionwe evaluate the performance of our 
method with respect to efficiency and to its incompleteness on several hundred 
examples originating from four different chemical and biological models. 


2 Finding Roots with Positive Coordinates 

Denote Ni = N \ {0}, and let d G Ni. For a G K, vectors x = (cci,... ,Xd) of 
either indeterminates or real numbers, and p = (pi, ... ,pd) G we use the 
notations ,..., a^'^) and We will, however, never 

consider a vector to the power of a number. Our notations are compatible with 
the standard scalar product as follows: 

(aP)'i = {aP^ ,..., = aP^i. 

Consider a multivariate integer polynomial 

/= coeff(/,p) • xP G Z[x], 

pGsupp(/) 

where coeff(/, p) 0 for p G supp(/), which is called the support of /. 

2.1 Finding a Point with Positive Valne 

The Newton polytope of / is the convex hull of supp(/). It forms a polyhedron 
in which we identify with its vertices, formally newton(/) C supp(/). The 
following lemma is a straightforward consequence of the convex hull property. 

Lemma 1. Let f = coefF(f,p) • xP + f' G Zfxl. Assume that p d newton!f). 
Then newton(/) = newton(/'). □ 
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For p e supp(/) we define sign(/, p) = sign(coeff(/, p)). We partition the 
support of / as follows: 


supp(/) 

supp+(/) 

supp"(/) 

supp°(/) 


supp+(/) U supp"(/) U supp°(/), 

{ p G supp(/) I sign(/, p) > 0 A p 7 ^ 0 }, 
{ p G supp(/) I sign(/, p) < 0 A p 7 ^ 0 }, 
supp(/) n { 0 }. 


Let supp+(/) = {pi,..., Pr}, supp (/) = {Pr+i, ■ ■ •, Ps}, and fix any order on 
supp(/). The basic LP matrix B{f) is composed as follows, where the last row 
is present if and only if supp°(/) yf 0: 


BU) 


B+{f) 


Pii 

Prd • • 

Pld -1 

■ • Prd 1 


= 

Pr+1,1 

■ • Pr-\-l,d 1 

B-{f) 


Ps,d 

Ps,d -1 

( 0 ,- 1 ) 


0 

0 -1 


Considering matrices concatenations of their rows, we write this also as B{f) = 
B^{f) o B~{f) o (0, —1)*. Whenever we write for a given matrix B G a 

product N* B, then we implicitly agree that 


N* 


-10 0 0 
0 10 0 

0 0 10 


^ ^mxm 


That is, the multiplication N*B replaces the elements of the first row of B with 
their additive inverses. Similarly, —1 = (—1,...,—1)^ is generally a column 
matrix of suitable length. In these terms, we are going to consider systems 

N* ■ B{f) ■ < —1, where x = (n, c) G 

which can be rewritten as follows: 

Pin — c > 1 

Pin — c < —1, iG{2,...,s}. 

Lemma 2. Let / G Z[x]. Let n G and let c G M. Then the following are 
equivalent: 

(i) The hyperplane iL(x) defined by nx = c strictly separates the point pi 
from supp(/) \ {pi}, and the normal vector n is pointing from H(x) in 
direction pi. In particular, pi G newton(/). 

(ii) There is 0 < A G M s.t. N* ■ B{f) ■ (An, Ac)^ < —1. 
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Proof. Assume that (i) holds. The orientation of n is chosen such that n-pi > c 
and n-pi < c for z S {2,..., s}. Define S = minig{i^..._s} |dist(pi, H)\ >0. Then 

Pi-n-c > (5|ln||, 

p,-n-c < -(^|ln||, te{2,...,s}, 

and we can choose A = ((5|jn||)“^. 

Vice versa, assume that (ii) holds. It follows that 

Pi • n > c+ 1/A, 

Pi • n < c-l/A, iG{2, ...,s}. 


Hence H{x) defined by nx = c is a hyperplane separating pi from supp(/) \ 
{Pi}, where the distance between H{x) and supp(/) is at least ||n|j/A > 0. 
Furthermore, n is oriented as required in (i). □ 


Lemma 3. Let 0 / G Z[x]. Then the following are equivalent: 

(i) There is (n, c) G s.t. N* ■ B{f) ■ (n, c)^ < —1. 

(ii) There is (n, c) G s.t. N* ■ B{f) ■ (n, c)^ < —1. 

(iii) There is n G c G Q s.t. N* ■ B{f) ■ (n,c)^ < —1. 


Proof. The existence of a real solution in (i) and a rational solution in (ii) 
coincide due to the Linear Tarski Principle: Ordered fields admit quantifier 
elimination for linear formulas m- Given a solution (ni,..., rid, c) G in 

(ii) , we can use the principal denominator m G Ni of ni, ..., rid to obtain a 
solution (mni,..., mud, me + m — 1) G Z'^ x Q in (iii). The implication from 

(iii) to (i) is trivial. □ 

Lemma 4. Let / G Z[x]\Z. Let (n, c) G such that A^*-i3(/)-(n, c)^ <—1. 
Then there is oq G N such that for all a G N with a > ag the following hold: 


(i) |coeff(/,pi) • 


> 


E 

z=2 


coeff(/,pj) • 


(ii) sign(/(a")) = sign(/,pi). 

Proof, (i) From / ^ Z it follows that pi 0. By Lemma we know npi > c 
and npi < c for z G {2,..., s}. It follows that there is 0 < <5 G K such that 


npi > c + S, (1) 

np* < c-5 zG{ 2, ...,s}. (2) 

We are going to show that oq = [max{2, {b ■ {k — 1)) ’’ }] is a suitable choice, 
where 

5 = |coeff(/,pi)r^ • max |coeff(/, Pi)|. 

iG{2,....s} 

For a > oo > 2 and for all z G {2,..., s}, the inequalities (Q and ([^ and 
monotony yield 

a"Pi > > a'*a"P’ > b ■ {k - 1) ■ a"P'. 
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function find-positive (/) 
data : / e Z[xi,. ..,Xd] 
result: p G or "failed" 

B+ -B+U) 

B- := B-{f) 
h := "infeasible" 

while h = "infeasible" and ^ \ ] do 
h := Ipsolve {B'^ o B~ o (0, —1)*) 
delete the first row from B~^ 

if h — "infeasible" then 
^ return "failed" 

(n, c) ;= h 
t := 2 

while /(t") < 0 do 
^ t:=2t 

return 

function Ipsolve (i?) 

data : a matrix B with d + 1 columns 
result: (n, c) G x Q or "infeasible" 

n := LP problem given by N*B and —1 
h := a solution (n,c) G of 11 or "infeasible" 

if h = "infeasible" then 
^ return "infeasible" 

m := principal denominator of the coordinates of n 
n := m • n 
c := me + m — 1 

return (n, c) 

Algorithm 1: Functions find-positive and Ipsolve 


Using the triangle inequality it follows that 

S S 

> |coeff(/,pi)r' • ^coefr(/,p,) • 

2—2 2—2 

which straightforwardly implies 


|coeff(/,pi) • a"Pi| > 


^coefr(/,p*) 


a 


npi 


2=2 

(ii) It follows from (i) that for a > ao the sign of the monomial coefT(/, pi) • 
^npi determines the sign of /(a”). Since a > 0, we obtain 

sign(/(a")) = sign(coeff(/,pi) • a"Pi) =sign(/,pi). □ 

After these preparations we can state our first subalgorithm as Algorithm 

Theorem 5 (Correctness of find-positive). Consider f G Z[x]. Then the 
following hold: 
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(i) The function find-positive terminates. 

(a) The function find-positive returns either "failed" or p G with 

/(P) > 0. 

Proof, (i) The termination of Ipsolve follows from the existence of terminat¬ 
ing algorithms for linear programming in line 15, including the Simplex algo¬ 
rithm [S], the ellipsoid method [IH], and the interior point method [17]. For 
the function find-positive itself, the number of iterations of the while-loop 
in line 3 is bounded by the number of rows of , which is in turn bounded by 
the finite cardinality of supp(/). The termination of the while-loop in line 11 
will be discussed with the correctness in (ii). 

(ii) To start with, the subroutine Ipsolve solves the LP problem II defined 
in line 14 and, in the feasible case, (n, c) in line 15 is a feasible point in 
The return value (n, c) G x Q in line 21 is a feasible point for II as well. Its 
construction in lines 18-20 follows the proof step from (ii) to (iii) in Lemma 
Next, the while-loop in line 3 has the following loop invariants. Consider 

n—1 

/(n) = f-Yl 


before the n-th iteration: 

(11) newton(/(„)) = newton(/), 

(1 2 ) !?(/(„)) =i?+)Oi3-o(0,-l)L 

Invariant (I 2 ) is easy to see. Consider (Ii). For n = 1 this is trivial. Before the 
n + 1-st iteration we know that h = "infeasible", which means that the LP 
problem given by 

7V*-((p„)oi?+^^)Oi3-o(0,-l)*) and -1 

was infeasible at the n-th iteration. According to Lemma [^ it follows that 
p>n ^ newton(/(„)). Using Lemma and the induction hypothesis we conclude 

newton(/(„+i)) = newton(/(„) - coeff(/, p„)xP’*) = newton(/(„)) = newton(/). 

The function find-positive has two possible exit points at lines 8 and 13 
corresponding to its two possible return values. Assume we are in line 13. We 
have to show that f" S with /(t“) > 0. The while-loop in line 3 has 

terminated after n iterations, and the if-condition in line 7 is false. In line 9 we 
know by (Ii), (I 2 ), and Lemma l^that the feasible regions for N* ■ o B o 
(0, —1)*) • V < —1 and N* ■ B{f) ■ v < —1 are identical. This allows us to use 
B'^ o B~ o (0, —1)* instead of B{f) for applying Lemmaj^to the original /, and 
our n S Z"^ has the property described there. In line 11, at the beginning of the 
k-th. iteration of the while-loop we have t = 2^. By Lemma we know that we 
will eventually have t> oq and thus sign(/(f")) = sign(/, p„) >0. □ 
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function find-zero (/) 
data : / e Z[xi,. ..,Xd] 
result: z G or "failed" 

V ■= /(I) 
if j/ = 0 then 
^ return 1 

if y > 0 then 

L / := -/ 

p:= find-positive(/) 
if q = "failed" then 
^ return "failed" 

z := construct-zero(/, p, 1) 
return z 
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function construct-zero (/, p, q) 
data : / G Z[xi,... ,a;d], p, q G 
result: z G or "failed" 

b := p + y • (q — p), where y is a new variable 
9 ■■= /(b) 

isolate r G ]0,1[ with g{r) = 0 
z := b(r) 

return z 


Algorithm 2: Functions find-zero and construct-zero 


2.2 Finding a Zero 

We have discussed how to heuristically find p G such that /(p) > 0 for 

our given / G Z[x]. On that basis Algorithm computes z G such that 

/(z) = 0, where Q denotes the algebraic closure of Q. 

Lemma 6 (Correctness of construct-zero). Consider / G Z[x], and let p, 
q G such that /(p)/(q) < 0. Then the following hold: 

(i) The function construct-zero terminates. 

(ii) The function construct-zero returns either "failed" or z G with 
/(z) = 0. If p, q G (Q'*")'^, then z G (Q''")'^. 

Proof, (i) The termination of construct-zero follows from the existence of 
terminating algorithms for univariate real root isolation including Sturm se¬ 
quences m and more efficient algorithms mm based on Vincent’s Theorem m- 
(ii) Since / is continuous and /(p)/(q) < 0, the intermediate value theorem 
guarantees the existence of z G pq with /(z) = 0. Formally, z G Q’’* is a solution 
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for X of the following nonlinear system with indeterminates xi, ..., Xd, y- 

/ = 0 (3) 

xi = pi+y{qi-pi) (4) 

Xd = Pd + y-{qd-Pd) (5) 

2 / > 0 ( 6 ) 

y < I- (7) 


In line 11, b is assigned the vector of the right hand sides of the d equations 

In line 12, these are plugged into the left hand side of equation ^ yielding 
a nonlinear univariate polynomial equation in y. Using any of the methods 
mentioned in (i), we obtain in line 13 a solution r € Q for y of that equation 
subject to the constraints That solution r is a real algebraic number 

in some suitable representation [23]. In line 14 we substitute r back into the 
equations Q“(j^ to finally obtain x = z G Q, also as a real algebraic number. 
Since z G pq, it follows from p, q G that also z G □ 

On the basis of Lemma the following theorem is straightforward. 

Theorem 7 (Correctness of find-zero). Let f G Z[x]. Then the function 
find-zero terminates and returns either "failed" or z G with /(z) = 0. 

□ 

When one is interested only in the existence of a zero of /, then one can, in 
the positive case, obviously skip construct-zero and exit from find-zero after 
line 8. Notice that, in addition, one can then also exit early from find-positive 
after line 8 in Algorithm 

2.3 An Illustrating Example 

Consider / = —2a;f + x\x 2 — 3a;f — X 2 + 2x2 ^ Z[xi,X 2 ]. We apply find-zero 
to find a point on the variety of /. Figure pictures the variety. We obtain 
/(1,1) = —3 < 0, and apply find-positive to /. 

Figure [^pictures the support of / and indicates the Newton polytope. We 
split into supp+(/) = {(2,1), (0, 2)}, supp-(/) = {(2, 0), (5,0), (0,3)}, and 
supp°(/) = 0, and we construct 
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Our first LP problem 











Figure 1: The variety of / = —2x1 + * 12^2 — 3xf — + 2x1 and the segment 

given by t G [0,2] of the moment curve corresponding to the normal 

vector (—3, —2) of the separating hyperplane in Figure]^ 



Figure 2: A subtropical view on / = —2xf + xlx 2 — 3xf — X 2 + 2x2 from Fig- 
ure[^ We see a hyperplane separating (0, 2) G newton(/) from supp(/)\{(0,2)} 
together with its normal vector n = (—3, —2). 


is infeasible, which confirms the observation in Figurej^that (2,1) ^ newton(/). 
Our next LP problem 
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is feasible with n = (-3,-2) and c = 5. Figure shows the corresponding 
hyperplane H given by —3xi — 2x2 + 5 = 0. It strictly separates (0, 2) G 
newton(/) from supp(/)\{(0, 2)}, and its normal vector n = (—3, —2) is oriented 
towards (0,2). We now know that > 0 for sufficiently large positive 

t. In fact, already 


/(2-3,2-2) = / 



1087 

16384' 
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The relevant part of the moment curve for t G [1,2] is pictured in 

Figure [2 Since both coordinates of n happen to be negative, the curve will for 
t —>■ oo not extend to infinity but converge to the origin. In particular, the curve 
will not leave the sign invariant region containing j). 

Finally, we call construct-zero with and (1,1) to solve the system 


—2x1 + ^ 1*2 — — X2 + 2x2 = 

Xi = 

X2 = 

y > 
y < 


0 


l + y (1 

0 

1 . 


I) 

i) 


Dropping a positive integer denominator, we obtain the univariate polynomial 

g = -16807J/5 - 12005?/4 - 934^^ _ 20778j/2 + 285?/ + 1087 

and an isolating interval y G ]0.2, 0.3[. Substitution of the real algebraic number 
( 5 , ]0.2, 0.3[) into the equations for Xi and X 2 yields an exact solution 

xi = ( 6861 ®-78a:3 + 584a;2 - 150a;-13, ]0.32,0.33[), 

X 2 = (l6807a;^- 12005a;‘‘ + 2026x3+ 9122x2-4609a:+ 323,]0.42,0.43[), 

where the intervals can, of course, be refined to arbitrary precision. Geometri¬ 
cally, our solving has intersected the variety with the line segment connecting 
the end points of our moment curve segment, which is also indicated in Figure 


2.4 Why Strictly Positive Coordinates? 

In the present section, we have focused on roots with strictly positive coor¬ 
dinates. This not only slightly simplifies the presentation. In fact, it is an 
important feature of our algorithm to be able to perform such a directed search. 

To start with, the research presented here was originally motivated by ques¬ 
tions on the stability of chemical and biological reaction networks, where the 
variables of the models typically are strictly positive. Our practical computa¬ 
tions in Section [5] are taken from those areas. For details on the theoretical 
background we refer the reader to [51 E51E51I551 15] . 

Furthermore, the concept of positive feasible points is well-known from linear 
programming. Techniques used there can be straightforwardly transfered to 
our situation: Consider / S Z[xi,..., x^j. For finding zeros (zi,, Zd) with 
sign( 2 :j) = Si G { — 1,1} consider /(siXi,..., s^x^), for zi G ]q;,oo[ consider 
/(xi + a, X 2 ,..., Xd), for zi S ] — 00 , /?[ consider /(—xi — /3, X 2 ,..., Xd), and for 
xi unbounded consider /(xi — x}, X 2 ,..., Xd) introducing a new variable x{. 


3 Finding Arbitrary Roots 

3.1 Using a Transformation 

Consider / G Z[xi,..., x^j. At the end of the previous section we have addressed 
a technique for turning a real feasibility test based on positive variables into a 
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general one. Using the observation that every real number is a difference of 
two positive real numbers, on introduces additional variables x'l, x'^ and 
transforms / into f{xi — x'l,... ,Xd — x'j). That transformation is ubiquitous in 
linear programming, if not explicitly then implicitly within the solvers. 

From an efficiency point of view our procedure is clearly dominated by the 
LP solving steps, where we have d+1 variables and |supp(/)| many constraints. 
Thinking in terms of state-of-the-art LP solvers [niiia and the Simplex method 
with the option of dualization m HO], the crucial complexity parameter is 
minjfi, |supp(/)|}. With our considered transformation the cardinality of the 
support increases exponentially in d in the worst case, but the number of vari¬ 
ables only doubles from d to 2d. 

Recall that our incomplete method relies on finding some p G newton (/) 
with coeff(/, p) > 0. We would like to point out that, doubling the dimension d 
with the transformation, the ratio |newton(/)| / |supp(/)| will in general increase 
for geometric reasons [S]. Furthermore the exponential increase of |supp(/)| 
increases the absolute number of candidates for a suitable p. On the other hand, 
the transformation does not add points to supp(/) but exchanges it entirely. It 
would require either comprehensive empirical studies or a thorough average-case 
analysis to make a precise statement about the quality of the transformation in 
terms of incompleteness. 

3.2 A Genuine Generalization 

We are now going to describe a generalization of the function find-positive 
in Algorithmic which searches for a suitable p not only in supp“*■(/) but also in 
a subset of supp“(/). Recall that in case of success, find-positive identifies 
the first row of B'^ as corresponding to the exponent vector of a monomial with 
positive coefficient that dominates / in the sense of Lemma|C Then it constructs 
a point p with large suitably balanced positive coordinates. The key idea for 
our generalization is the following: If the coefficient of an otherwise suitable 
monomial is negative but there is at least one odd exponent in the exponent 
vector, then we can correct the “wrong” sign of the coefficient by replacing the 
respective coordinate in the constructed point p with its additive inverse. 

For p = (pi,... ,pd) € supp(/) define the minimal odd coordinate 

moc(p) = min{i G {!,..., d} | 2\pi}, 

where min 0 = oo. We use the minimal odd coordinate to partition supp“(/) = 
supp“’“(/) Usupp®“, where 

supp“'"(/) = { p G supp" I moc(p) < oo }, 
supp'*"(/) = { p G supp" I moc(p) = oo }. 

The elements of supp'“'"(/) are called weakly negative. They have at least one 
odd coordinate. The elements of supp'*"(/) are called strongly negative. They 
have exclusively even coordinates. We furthermore define and B^~{f) 

corresponding to supp'“'"(/) and supp®"(/), respectively, and we obtain 

B{f) = R+(/) o R-"(/) o B^-{f) o (0, -1)U 

Consider a matrix B obtained from B(f) by deleting rows. Then we define 
moc(R) = moc(Rii,..., Bid), i-e., the minimal odd coordinate of p G supp(/) 
corresponding to the first row of B. 


11 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 

22 


function find-positive-general(/) 
data : / e Z[xi,. ..,Xd] 
result: p G or "failed" 

B-.= B+{f)oB^-{f) 

B^- :=B*-(/) 

/r := oo 

h := "infeasible" 

while h = "infeasible" and i? 7 ^ [ ] do 
if the first row of B is in B'^~{f) then 
^ ^ moc(i3) 

h := Ipsolve {B'^ oi?o( 0 ,— 1 )*) 
delete the first row from B 

if ft, = "infeasible" then 
^ return "failed" 

(m ,... ,nd,c) := h 
t := 2 

:= (r%...,r0 

if /r < 00 then 

^ til ■= —tfi 

while /(ti, ...,td) <0 do 
t := 2t 

if < 00 then 

^ til '■= —tfi 

_ return {ti,... ,td) 

Algorithm 3: Function find-positive-general 


After these preparations we can state our function find-positive-general 
in Algorithm]^ A corresponding function find-zero-general is obtained by 
replacing in find-zero in Algorithm the call to find-positive with a call 
to find-positive-general. Everything else remains unchanged. 

For showing the correctness of find-positive-general we are going to use 
the following variant of Lemma ^ 

Lemma 8. Let / G Z[x]\Z. Let (n, c) G such that A*-i3(/)-(n, c)^ <—1, 
and let fi = moc{B{f)) < 00 . Then there is oq G N such that for all a G N 
with a > Qq the following holds: Define t = (ti,..., G with tj = a”^ for 
j G {1,... ,d} \ {n} and Oii = —t"''. Then 

sign(/(t)) = -sign(/,pi). 

Proof. We have tj = a"j > 0 for j G {1,..., d} \ {/r}, = —a"'" < 0, and 

is odd by definition of the minimal odd coordinate. It follows that 

0 <a"Pi=-tPL ( 8 ) 

For i G {2,...,s} we have at least |a"P‘| = |tP’|. This allows us to conclude 
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from Lemma |^(i) that 


|coefr(/,pi) • tPi| > 


^coefr(/,p*) -tP’ 


z=2 

Hence coefT(/,pi) • tP^ determines the sign of /(t). Using the inequality in ^ 
we obtain 


sign(/(t)) = sign(coeff(/,pi) -tPi) = -sign(/,pi). □ 

Theorem 9 (Correctness of find-positive-general). Consider f G Z[x]. 
Then the following hold: 

(i) If the function find-positive in Algorithm^ does not fail on f, then 
find-positive-general(/) = find-positive(/). 

(a) T/ie/wnch'on find-positive-general terminates. 

(Hi) The function find-positive-general returns either "failed" or p G 
with /(p) > 0. 

Proof, (i) The function find-positive-general operates on B~^{f) oB'"~ {f) o 
B^~if) o (0,-1)* while find-positive operates on -B+(/) o B~ o (0,-1)* 
so that there is possibly a different order of rows lying below i?+(/). How¬ 
ever, when find-positive does not fail, then the same feasible solution is 
found in both functions before touching anything outside B^, and in line 10 of 
find-positive-general we have exited the while-loop with /i = oo. It follows 
that the if-conditions in lines 15 and 20 of find-positive-general are always 
false, and the rest of the code after the while-loop is computationally equivalent 
to the corresponding part of find-positive except for an expanded notation. 

Accordingly, a proof of parts (ii) and (iii) can be straightforwardly de¬ 
rived from the proof of (i) and (ii) of Theorem respectively: If the function 
find-positive in Algorithm 1 does not fail on /, then there is nothing else to 
do. Otherwise we always reach lines 15 and 20 with fj, < oo, replace t^ with its 
additive inverse, and apply Lemma instead of Lemma [^(ii), where we know 
that that sign(/, p„) <0. □ 

4 Practical Issues 

In this section we would like to discuss issues and share experiences related to 
a practical implementation of our method. 

One major benefit of our approach is the reduction of an algebraic problem to 
linear programming (LP). Linear programming is a field with more than 50 years 
of active algorithmic research, strongly driven by practical applicability and 
aiming at robust implementations. Our own implementation combines the the 
Codemist Standard Lisp (CSL)-based version of the computer algebra system 
Reduce naniiis] with the Gurobi Optimizer |13j . Technically, CSL provides 
a foreign function interface that allows us to dynamically load the Gurobi C- 
library at runtime and call its functions from within symbolic mode Reduce 
functions. Gurobi uses the Simplex algorithm. So far we have got no experience 
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with the use of implementations of polynomial methods for LP, like the interior 
point method m- 

Gurobi uses floating point arithmetic with a limited precision. We want to 
adress some issues related to this, which we consider of general inteterst, because 
that floating point approach is typical for Simplex-based LP software. 

On the one hand, LP solvers are quite good at controlling numerical stability. 
With our comprehensive computations we have never encountered any problems 
with false results due to LP rounding errors. On the other hand, in line 15 of 
Algorithmj^we obtain ni, ..., n^, c as floats with small rounding errors. These 
rounding errors do not affect correctness but cause a subtle problem: Convert¬ 
ing Til, ..., Ud into fractions, the GCD of their denominators will typically be 
1 so that the principal denominator m computed in line 18 becomes the very 
large product of those denominators. Consequently, in line 19 we obtain our 
final n with very large relatively prime integer coordinates. This, in turn, ren¬ 
ders infeasible the exponentiation of increasing powers of 2 with those integer 
coordinates and substitution of the result into / in line 11 of Algorithm or in 
lines 14, 17, and 19 of Algorithm]^ There are two principle ways out, which we 
call the pure LP approach and the MIP approach, respectively. Of course, the 
single design decisions made with these approaches can be recombined to yield 
further, mixed, approaches. 

The Pure LP Approach The pure LP approach tries to get along with the 
delivered floats. Specifically, lines 18-20 in Algorithm are skipped, and a 
floating point vector is returned. The while-loops in line 11 of Algorithm and 
line 17 of Algorithm remain correct with floating point exponents n. Later, in 
lines 11-12 of Algorithmit is important to convert to rationale. In particular 
the substitution of floats into a high-degree polynomial / in line 12 could cause 
considerable numerical instabilities. Subsequent root isolation to floating point 
precision in line 13 and back-substitution of the obtained floats in line 14 worked 
well with all our computations. 

The MIP Approach MIP stands for mixed integer (linear) programming. 
Our Lemma allows us to declare ni, ..., as integers to the LP solver right 
away, while c remains real. As MIP is NP-hard m. the MIP approach is con¬ 
siderably harder than the pure LP approach in terms of theoretical complexity. 
In practice there are several advanced algorithms for Simplex-based MIP solv¬ 
ing, which rely in some way on considering an LP relaxation, i.e., considering 
integer variables as real variables, and, in the feasible case, trying to construct 
a mixed real integer feasible point on the basis of an LP solution. The Gurobi 
solver specifically uses advanced cutting plane m methods for that construc¬ 
tion. For the largest problems discussed with our practical computations in 
Section below, we have observed factor of about 3 for MILP solving compared 
to LP solving. 

There is an interesting optimization with the MIP approach: Since in out 
situation MIP feasibility is equivalent to LP feasibility by Lemma one can 
generally first check the latter in lines 14-15 of Algorithm and in the feasible 
case rerun for the corresponding MIP problem. Using this strategy, there is 
always at most one MIP solving step per root finding problem. Furthermore, one 
runs MIP solving only on feasible instances. This excludes the really problematic 


14 


13 {yi,...,yd) = ,..., 2/”'^) for a new variable y 

14 if /I < oo then 

15 L 

16 /i :=/(yi,...,?/d) G Z(2/) 

17 t := 2 

18 while fi(t) < 0 do 

19 L ^ 

20 return 1 /d(t)) 

Algorithm 4: Code to replace lines 13-22 in Algorithm 


cases, which are LP feasible but not MIP feasible problems. 

In rare cases one obtains integer solutions which are so large that they render 
exponentiation and substitution in line 11 of Algorithm or in lines 14, 17, and 
19 of Algorithm 1^ infeasible. One can impose a suitable bound on the absolute 
values of the solutions, and in case of exceeding that bound treat the problem 
as infeasible, and proceed to the next candidate. 

Another noteworthy optimization is the symbolic precomputation of a uni¬ 
variate rational function for the while-loop in line 17 of Algorithm See Algo¬ 
rithm for details. A corresponding simpler variant, of course, works also for 
lines 10-12 in Algorithmic 

For root isolation in line 13 of Algorithm |C we use the Vincent-Collins- 
Akritas method [1] . We obtain a real algebraic number encoded by a univariate 
defining polynomial and an open isolating interval, which is back-substituted in 
line 14, yielding a vector of such real algebraic numbers as the final solution z. 

5 Some Practical Computations 

We consider input polynomials originating from 4 different chemical and bio¬ 
logical models. This yields 929 instances altogether. For all of these instances 
we are checking for zeros with strictly positive coordinates. It turns out that 
for 640 of the instances we find = [ ] in line 1 of Algorithm [C which tells 
us that the corresponding polynomial is positive definite (on the interior of the 
first hyperoctant). Running our method on the 289 remaining instances, it fails 
in only 7.3 percent of the cases. Table shows detailed information for the 
single models. It also shows size (number of monomials), dimension (number 
of variables), and the largest degree of an occurring variable for the respective 
largest instance. It furthermore shows the maximal computation time for a sin¬ 
gle instance and the sum of computation timesQ All computations have been 
carried out on a 2.8 GHz Xeon E5-4640 with the MIP approach, yielding exact 
algebraic number solutions: 

Notice that for our particular application the detection of definiteness by 
our implementation establishes a perfect result. From that point of view, one 
could argue that our method fails in only 3 percent of the cases. 

^All input and log files are available at http://research-data.redlog.eu/arXiv/2015/ 
subtropical/ 
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METH 

OMBO 

MBO 

MAPK 

Total 

number of instances 

7 

496 

405 

21 

929 

number of definite instances 

3 

338 

283 

16 

640 

number of remaining instances 

4 

158 

122 

5 

289 

found zero in 

4 

144 

107 

5 

260 

failed on 

0 

14 

15 

0 

29 

failed on (% of remaining) 

0 

8.9 

12.3 

0 

7.3 

size of largest instance 

347 

9787 

9706 

863438 

863438 

dimension of largest instance 

7 

7 

7 

10 

10 

degree of largest instance 

6 

10 

9 

12 

12 

maximal time (s) 

0.16 

4.68 

10.00 

15.87 

15.87 

total time (s) 

0.21 

199.91 

162.88 

15.92 

379.92 


Table 1: Statistics for our practical computations 
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